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1.0 Introduction and Summary 


1.1 Overview 


Geomagnetic field models are important for a number of uses, includ- 
ing investigation of processes deep within the earth's interior, trend 
removal for anomaly studies, and the investigation of the correlation 
between geomagnetic and geodetic anomalies. (Field models are also used 
for navigation and survey purposes and charged particle trajectory calcu- 
lations.) The typical method of modeling the main field and the secular 
variation is by least squares analysis using spherical harmonic models. 

An extensive review of spherical harmonic methods and main field and sec- 
ular variation models is given by Barraclough [1975, 1976']. 


Most users of geomagnetic field models require accurate estimates of 
the current field. This necessitates prediction beyond the data domain 
of a field model and, inevitably, involves extrapolation inaccuracies 
caused by uncertainties in modeling changes in secular variation 
patterns. Due to present lack of knowledge of the underlying physical 
processes involved, these patterns can typically be accurately modeled by 
linear variations only over periods on the order of a few years. 


Recursive estimation theory provides a means of combining models 
obtained by conventional batch least squares over deterministic periods 
into optimal estimates of the field model parameters at any particular 
time. In particular, this technique should provide more accurate field 
model prediction capability. This is done by using statistical informa- 
tion about the temporal variation in the field model to weight the indi- 
vidual least squares solutions for a particular time interval (here 
referred to as a batch estimate). Obviously, the most recent data is 
weighted most heavily and past data is given the least weight. However, 
the data from all of the conventional solutions over the deterministic 
periods is included in the solution so that the prediction errors should 
be significantly less than those for a single batch estimate. In 


' ^1 
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particular, the MAGSAT mission, with accurate vector and scalar data, 
will provide for a highly accurate field model for epoch 1980, and the 
recursive procedure will allow for improved prediction capability from 
the MAGSAT model by optimally incorporating past information. 

Other advantages of recursive estimation include the ability to 
easily update the field model as new data becomes available, improve the 
field estimate for times in the past and also to generate an estimate of 
the error in the field. 

1.2 Summary 


The results of a preliminary study to determine the feasibility of 
using Kalman filter techniques for geomagnetic field modeling are given 
in tins report. Specifically, five separate field models were computed 
using observatory annual means, satellite, survey and airborne data for 
the years 1950 to 1976. Each of the individual field models used approx- 
imately five years of data. Then these five models were combined using a 
recursive information filter (a Kalman filter written in terms of infor- 
mation matrices rather than covariance matrices.) The resulting estimate 
of the geomagnetic field and its secular variation was propogated four 
years past the data to the time of the MAGSAT data. The accuracy with 
‘which this field model matched the MAGSAT data was evaluated by compari- 
sons with predictions from other pre-MAGSAT field models. The field 
estimate obtained by recursive estimation was found to be superior to all 
other models. 

It must be emphasized that this study was not intended to be the 
final word on field modeling. Because of the preliminary nature of the 
study, several "short cuts" were taken which make the field estimate sub- 
optimal. In particular, Kalman filtering can only be optimal when the 
statistics of the unmodeled field dynamics (among, other quantities) are 
known. In this study, there statistics were obtained by an "eyeball 
estimate" using plots of the estimated field model coefficients and were 
refined using the likelihood function computed by the Kalman filter. 
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However, these estimates are crude at best and should really be obtained 
using a more sophisticated technique such as maximum likelihood estima- 
tion. (Obviously it would be preferable to use knowledge of the core 
dynamics but this is not well known at the present time). 

Another "short cut" involves the handling of observatory "biases" 
and scaling of the information matrices to account for aliasing and data 
noise. This is discussed in detail in later sections. ' 

Considering the limitations of this study, we believe that the 
results are encouraging and warrant further analyses to optimally apply 
recursive estimation for geomagnetic field modeling. 
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2.0 Mathematical Description 

2.1 Deterministic Geomagnetic Models 


While several organizations are currently involved in core field 
modeling, the basic techniques used date back to Gauss, with virtually no 
utilization of modern estimation methods. Typically, the magnetic field 
is parameterized as .. scalar potential expanded in spherical harmonics 


” ^internal ^ ^external 


j;(« 

k /r\^ n 

Mi ^ 

n=l'®^ m=( 


I (g|^cosm'|)+h|!|sinm<l>) P^(0) 
m=0 


(Y|!|cosmi{i+tS^sinm<J>) p||I|(Q) 


where 

a = mean radius of the earth 
n* = maximum degree of the expansion 

r = geocentric distance 

g, h = quasi-normalized Schmidt coefficients 

ifi = east longitude 

9 = colatitude 


P^(0) = associated Legendre Functions (Schmidt quasi-normal ized) . 
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Here V is the potential at a point in a domain free from sources of 
magnetic field, and the field vector B is given by the gradient of this 
potential 


B = -VV . 

The coefficients g and h are typically represented deterministically as 
first or second order Taylor expansions in time about some epoch t^, with 
a conventional least squares batch processing method (outlined in Section 
2.3) used for their recovery. 

While data from repeat stations, marine surveys, aeroinagnetic 
surveys, magnetic observatories and satellites (principally the POGO 
series and currently MAGSAT) contribute to the determination of the geo- 
magnetic field over a time interval about the epoch, annual means from a 
world-wide network of magnetic observatories (see Figure 2.1) represent 
the most useful data set for determining the secular variation of the 
internal field. The difficulty in modeling the temporal variation of the 
field over more than a short time interval is illustrated in Figure 2.2 
for the Abisko observatory. There the solid lines represent a 
quadratic polynomial model for g^ and h^ over the interval 1900-1965. 

For prediction of the magnetic field, the selection of the temporal 
representation for a given data span is crucial, and the danger of 
extrapolating beyond the data is clear. The incorporation of annual 
means data into a main field model for a particular epoch suffers from 
the fact that the magnetic field measured at the observatory is the 
vector sum of the main field and a contribution due to local crustal 
magnetization, 


^ "* ^internal ^ ^crustal 


where may change appreciably over the distance of a few 

kilometers. While varies with time, however, remains 

constant. Thus, models of secular variation based only on time 
derivatives of annual means observations 









1 


Business. 4M) TECHmunmi tL Systems, /yvr. 







'rt 



'i 


^ “ ^internal 

are not influenced by the local anomalous fields, but models of the main 
field are aliased. 

In this analyses, the main field and its secular variation were 
estimated simultaneously, using both annual means and other data types 
(satellite, airborne, survey). The annual means data is accommodated 
by solving for the local observatory bias, , at each observa- 

tory. This allows the data to properly distribute its influence among 
the secular variation and constant parameters of the model in a least 
squares sense. The local biases, which are estimated along with the 
field model , provide some physical measure of the local anomaly field. 

2.2 Modeling The Temporal Variation of the Field 


In the previous section, it was stated that the temporal variation 
of the core geomagnetic field has usually been modeled deterministically 
as a first or second order Taylor series expansion. This deterministic 
modeling is a necessity when conventional least squares estimation is 
employed. However, when Kalman filtering is used, the process can be 
modeled as a stochastic one. Thus the question arises as to whether some 
form other than a truncated Taylor series would be more appropriate for 
modeling the temporal variation. For example, a first or second order 
Markov process may be suitable. However, Markov process models are 
mainly useful when the data span is a significant fraction of the time 
constants of the process and when the process is zero mean. If the data 
span is much shorter than the time constants of the process, the output 
can be approximated quite well as the integration of initial conditions 
and white noise inputs for a system which does not use feedback (although 
the order of the system may be di f ferent than that of the Markov 
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process). In other words, knowledge of the time constants of a Markov 
process does not significantly improve the modeling of that process when 
the data span is much shorter than the time constants. 

For the case of geomagnetic core field modeling, the shortest time 
constants appear to be 50 to 100 years (it is difficult to determine them 
accurately when only 80 years of data are available). Since we are only 
using 26 years of data (1950-1976), the use of a Markov model is not 
justified (even if it was, we do not know the time constants). 
Furthermore, many of the low degree coefficients in the spherical 
harmonic expansion are definitely not zero mean (e.g. 9 jq)* Thus the 
"bias" in coefficients would have to be estimated separately from the 
Markov process, adding unnecessary complexity to the model. 

For these reasons, Markov models were not used in our analyses. 
Rather, the stochastic contribution of the temporal variation of the 
spherical harmonic coefficients was modeled as one of two forms: 

low degree: c’. (t) = w^(t) 

high degree: c. (t) = w^(t) 

where: 


c^. represents a coefficient of the spherical harmonic expansion 
w^ is random noise. 

Thus the expressions for c as function of time are: 

low degree: c.(t) = + c.(tQ)(t-tQ) + c. (tQ)(t-tQ)^/2 

+ q^(t-tQ) 

high degree: c.(t) = c.(tQ) + c-(tQ)(t-tQ) + q^-(t-tQ) 
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where q^(t-to) represents a random disturbing input (i.e. weighted 
integral of w(t) (See Section 2.3.3). 

The dividing line used in deciding between the linear and quadratic 
models was set so that the differences between the data and the model 
appeared to be random. All coefficients of degree 7 or less used the 
quadratic model while the linear model was used for higher degree 
coefficients. 

2.3 Estimation Procedure 


Data available for field modeling is comprised of survey and 
observatory annual means values (primarily D, 1, H, 2, F), oceanographic 
data (F), aeroinagnetic data (F), and satellite data (F). Witli MAGSAT , 
there will be available satellite vector component data as well. For 
this preliminary study, only data from 1950 to 1976 was used. Pi'ior to 
1950, the data coverage Is sparse. No data after 1976 was used so that 
the estimated field could be predicted to 1980 and compared with the 
MAGSAT data. During the period 1950 to 1976, approximately 186,000 
observatory, ground survey, aeromagnetic and satellite observations are 
available. Oceanographic and repeat survey data were not used in this 
preliminary analysis because of its questionable quality and difficulty 
in using it properly. 


The following sections describe the procedure used to compute 
"mini -batch" field models, the recursive estimation procedure and the 
calculation of the state noise matrix. 


t 2.3.1 Mini -Batch Field Models 

I 

S* 

I The recursive estimatir;n procedure used here requires the 

s 

I calculation of conventional epoch field models on a "mini-batch" interval 

I of approximately 5 years of data. The computer software used for this 

I 
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purpose is that originally developed at Goddard Space Flight Center by 
Cain, et al , 1967 and recently modified and enhanced by Business ana 
Technological Systems, Inc., (Estes, 1979 and 1980). 

For the purpose of clarifying equations to follow, we shall 
summarize the equations of the conventional batch solution in matrix 
notation. Values of the g and h coefficients at any point in time are 
assumed to be in terms of the values at epoch (tg). Since we are using 
only 5 years of data in each “mini-batch," the field can be adequately 
represented (with negligible aliasing) using only constant and linear 
terms, ie. 

g(t) = g(tg) + g(tQ)(t-tQ) 
h(t) = h(tQ) + h(tQ)(t-tQ) 

The maximum degree to which the constant and linear terms in the 
expansion was carried varied for the different mini-batches. When 
satellite data was available (after 1960), the constant and linear terms 
were carried to degree 13. However, prior to 1960, the available data 
would not support a solution with so many unknowns. Thus the cutoff was 
lower. iDetails for the individual mini-batches are given in 
section 3.1). 


For purposes of clarity in the discussion to follow, we will assume 
that secular velocity terms are present through n* = 13. The 
relationship of the parameters can be written in matrix form by defining 
a state vector which contains the values and derivatives of the 
coefficients. Let 


.. _ ,„1,0 -1,0 „1.1 *1.1 .1,1 J3,13 :13,13 .1 ^1 .1 

x=Lg’ g g g h’ h’ ...h h ’ b b b ...b J 




where b^ , b^, b^ are "biases" in the x, y and z components of the 

A V fc 


magnetic field at observatory i. {This accounts for localized fields) 
Then 


where 


x(t) = *(t-tQ) x(tQ) 
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Let the set of measurements during the mini -batch period be denoted 
by vector y of length in (ni may be on the order of tens of thousands). 
Then y can be written as a nonlinear function of x plus random measure- 
ment noise (y = f(x) + v). The perturbations in y (Sy) due to 
perturbations in x (5x) can be written as 

6y = M 5x 

where 

3y 

M = ^ (dimension mxn) . 

Then the weighted, least squares, differential correction estimate of 
is given by 


XI 
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where: 
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= the estimate of x at iteration K 


W = the covariance matrix of the random measurement errors. If 
y = f{x) + V where v is a zero-mean, random error vector, then 


W = E(vv ) (diagonal ). 


Denote the converged output of mini -batch number i as x^^^ . Then the 
associated covariance matrix for the error is given by 




This matrix is produced as a by-product of the least squares calculation. 


After convergence, those partitions of the Information matrix 
= P^l) and the state estimate Xj^. corresponding to the field 
model (spherical harmonic coefficients and derivatives) are written on 
tape for later use. However, this partition of the information matrix 
must be equal to the inverse of the corresponding partition of the 
covariance matrix. If the original information matrix and vector are 
partitioned as follows: 




, M'''w"^ay 


where the first partition corresponds to the spherical harmonic 
coefficients and the second partition corresponds to the observatory 
biases, then 
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-1 T -1 
-A3^A2D ^ 

-1 T 

where D = ^2 partition of the information matrix written 

on the tape. 

It was sometimes necessary to modify the epoch of the mini-batch 
after the infomiation had already been stored. This is done using the 
following transformations 

D(t) = <J.^(tQ-t) D(tQ) MtQ-t) 

X(t) = <|>(t-tQ) X(tg) 

2.3.2 Optimal Filtering of Mini-Batch Field Models 

The output from the mini -batches must be combined to produce an 
optimal estimate of the spherical harmonic coefficients at the epoch of 
the last mini-batch (1972.5). When properly done, this estimate should 
yield good predictions of the geomagnetic field for several years into 
the future. The key to combining the results of the mini -batches is to 
correctly weight their output. This can be done via a Kalman filter 
which includes the effects of process (state) noise (Gelb, 1974). 

Recall that the temporal model for the spherical harmonic 
coefficients was assumed to be 

x(t) = *(t-tQ) x(tp) 

where the definition of has been expanded (because of the longer time 
span) to include quadratic terms in addition to the constant and linear 
terms for coefficients of degree 7 or less. The partitions of 
corresponding to quadratic terms are 


-D’^AgAg^ 

A 3 + A 3 ^ A 2 D ^ 3 ^ 
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Although this model is quite adequate when the time span is short, it 
degrades rapidly as the time span grows. When combining the estimates 
from different mini-batches, some allowance must be made for the errors 
in the secular model (e.g., third and higher order derivatives). In 
other words the secular model can be better written as 

x(t) = f(t-tQ) x(tQ) + g(t-tQ) 

where g is assumed to be an unknown, random disturbance vector. 

We now change the notation so that time is not explicitly shown: 

-N+l = 5l ^ 3l • 

Here, the subscript refers to the epoch of each mini-batch. Since each 
mini -batch will have a different epoch {at the middle of the data span), 
an equation of this type must be used to relate the coefficients (state 
vector) for different mini-batches. However, since q^. is unknown 
(to us), our best estimate of x.^^ based upon the measurements included 
in mini -batch i is 

where means the estimate of x^. at mini-batch i based upon measure- 
ments up to mini-batch j. is called the a priori estimate at 

time i+1 and x.^^. is called the a posteriori estimate at time i. 


The a priori error covariance matrix for x... 


is easily shewn to 
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where Q^-=E(g^. gj). 


For the moment, we will ignore how is obtained. 


After mini -batch i+1 is processed, we can obtain an a posteriori 
estimate for which optimally combines the results of the new mini- 
batch and the results of previous mini -batches. 




where T is a matrix consisting of zeroes and ones which shifts the 
elements of to match those of x^^^ (x^ does not contain quadratic 

time terms). This is called an information filter (rather than a Kalman 
filter) because the result is obtained by combining information 
matrices. For this particular problem, it is more efficient than a 
Kalman filter. At each step, the inversion of two, positive-definite 
synsnetric matrices is required. This is best done by Cholesky 
factorization which takes a total of n^/2 operations. For the dimension 
of the problem discussed previously, this requires a relatively small 
amount of computer time. 


There are two aspects of this procedure which deserve further 
discussion. First, it was not necessary that all of the coefficient sets 
be the same size for all mini -batches. For example, it was not desirable 
(because of the sparseness of data) to solve for all high degree 
coefficients and derivatives for periods prior to 1960. This simply 
means that the partitions of the information matrix (M^W“^M) correspond- 
ing to these coefficients contained zeroes. Secondly, notice that new 
data can be easily incorporated into the solution without rerunning the 
entire procedure. As additional satellite data is obtained, it can be 
used to improve the field model. 
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The final output from this procedure is a state vector (coefficients) 
and error covariance matrix which are evaluated at the epoch of the last 
data span. These can then be used to predict the coefficients and error 
covariance for the next few years. This estimate should be better than 
the estimate obtained from the last mini -batch because it uses more data 
going back further in time. However, notice that old data is given less 
weight in the final solution than new data. The inclusion of Q tells the 
filter that the state vector is subject to random changes and, thus, it 
gives the most weight to recent data. 

The ability to estimate the error in the predicted coefficients is 
an important output of this procedure. Although all batch estimation 
procedures produce an error covariance matrix as a by-product, the use- 
fulness of this covariance has not been widely appreciated. Since the 
batch procedures do not include a Q term, their error covariance matrices 
cannot be used to accurately predict the field errors for times outside 
of the data span. However, the procedure proposed here can do this. By 
a simple transformation on the covariance matrix, it is also possible to 
compute the error in the estimated magnetic field at various points on 
the earth's surface. This could be used to produce a contour map of the 
earth showing the accuracy of the estimated field. 

The accuracy of the proposed estimation procedure depends upon the 
accuracy of the Q which is used. Past experience has shown that, for 
most problems where Q is small, the final estimate is not very sensitive 
to Q. If Q is known within an order of magnitude, the output will 
probably be satisfactory. 

The next section derives the mathematical equations used to compute 
Q while section .2.3 shows how the likelihood function can be computed in 
the Kalman filter and used to help refine the estimate of Q. 


The state noise covariance matrix, Q, is defined as the covariance 
matrix of the unmodeled dynamics in a discrete system. can be calcu- 
lated from the continuous system if the spectral density of the random 
input is known. In other words, if the continuous system is defined as 

x(t) = Fx(t) + u(t) 

T 

where u(t) is white noise, then can be computed from E(u(t)u (t)). 
Using this equation for the continuous system, the corresponding equation 
for the discrete system is: 

x(t) = ‘f(t-tQ)x(tQ) + <5(t-X)u(X)dA 

where the second term is equal to q^ as described in section 2.2.2. Thus 

t . t . 

q. = E «(t.-x)u(x)dx / i u'’'{x)$T(t.-x)dx 

^ S'-i ^ ^i-r ^ 


” t. T T 

= E / ' ♦(t.-X}u(X)u‘{X}4'(t.-X)dX 

_ ^ 

t . -r T 

/. *(t.-X)E[u(X)u'(X)] x)dX . 

^i-1 1 - - 1 

T 

It is assumed that E[u(X)u (X)] is a constant (Q^). For the polynomial 
representation of the spherical harmonic coefficients, the 4 and 
matrices take one of two forms: 
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quadratic terms: 
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linear terms: 


«(At) = 


0 0 

0 Q. 


When these expressions are substituted in the expression for Q- and 

• 1 

integrated, the following results are obtained: 


quadratic terms: 


^i ^s2 


1^5 

At^ 

At^ 

5 

2 

“T' 

At^ 

4At^ 

At' 

2 

3 


A4.3 

[At 

At^ 

At 


ORIGINAL PAGS IS 
OF POOR QUALITY 


linear terms: 


^i " ^sl 2 

1 Si 


We have implicitly assumed that the process noises for coefficients 
of different degree and order are uncorrelated. This seems to be a 
reasonable assumption considering the orthogonality of the spherical 
harmonic coefficients. 
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2.4 Computation of the Likelihood Function qP pooR QUALITY 

In any estimation problem where the measurement and dynamic models 
are not known with certainty, it is generally useful to compute the log 
likelihood function within the Kalman filter. This likelihood function 
can then be used to aid in defining the model. When filter runs using 
different model parameters are compared, the run which has the highest 
likelihood function can usually be selected as having the "best" model. 

In an ordinary Kalman filter, the extra burden required to compute 
the log likelihood is trivial since the required quantities are already 
computed by the filter. However, in our information form of the Kalman 
filter, the computational burden is quite significant. The total 
procedure for the information filter and computation of the likelihood is 
given below: 


Step 

Filter 

Computation of 
log likelihood 

1 

Read t^,(M^W“^M)^ and ^ 

' 

2‘ 

SVi-l " * ^i-l/i-1 



'"i/i-l = * ^-l/i-l * ^^i 


3 


’'i " Vi/i-1 ■ ^b,i 

4 


L = L+xlc"^xT+ln(det(C)) 

5 

Pi/i = 
’‘i/1 " 
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where: and Tg are transformation matrices consisting of zeroes and 

ones which change the dimension of the filter state and the 

mini-batch state to match. is used to increase the dimension 

of X. to match x and is used to decrease the dimension of x 
° . 2 

to match Xj^. In the computer program, this transformation is 
performed by shifting operations (no matrix multiplications are 
involved), 

L is equal to -2 log likelihood (plus a constant bias). 

T -1 

The inversion of (M W M). in step 3 and the inversion of C in step 
4 are the most costly operations in computing the log likelihood. In 
fact, these operations almost double the running time of the program. 
However, the information obtained from the likelihood function is so 
valuable that the extra burden is well worth the price. 

Since the combined operation of an information filter and some 
computations of the Kalman filter (steps 3 and 4) appears to be somewhat 
redundant, the question arises as to whether everything could be done 
more efficiently using just a Kalman filter or an information filter. 

The answer is no. In fact, if a Kalman filter were used to compute 
everything, the running time of the program would double. (This 
statement is not generally true for most filtering applications). 

For more information on the properties of the log likelihood 
function, see Edwards, 1972 and Gupta, 1974. 
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3.0 RESULTS 

3.1 Pre-MAGSAT Model PMAG (7/80) 

The models used as a basis of comparison for prediction accuracy 
include the AWC 75 (Peddie & Fabiano), the WC 80 (Barker & Barraclough) 
and the Goddard Space Flight Center Model, PMAG (7/80) (Estes). Because 
the PMAG (7/80) model is not fully documented in the literature and 
because the techniques utilized in the development are also used in 
deriving the mini -batch models of Section 3.2, a brief description of the 
modeling method will be presented in this section. 

The basic core field data set consisted of POGO satellite scalar 
data, selected observatory vector data, and selected repeat station data 
and marine survey data over the years 1960-1977. Magnetic potential 
spherical harmonic coefficients were simultaneously least squares fit to 
the data set to mathematical degree and order 13 for the constant terms, 
13 for the first time derivative terms, 6 for the second derivative 
terms, and 4 for the third derivative terms. While the satellite data 
provided a good uniform spatial distribution, vector component data was 
desired to reduce the Backus effect and observatory annual means data was 
needed to provide accurate secular variation information. The surface 
data, however, is the result of contributions from both the core field, 
B(t) , and the anomalous crustal field, , which can contribute 
several' hundred nT to the measurement but is assumed not to vary with 
time. To properly account for the crustal influence in the data, differ- 
ent approaches were taken for satellite and for surface observatory, 
marine, and repeat station data. 

The basic 49,000 measurement POGO data set used in developing the 
POGO 2/72 model (Langel, et al ) was extended in time to include 
approximately 25,000 additional OGO-6 quiet measurements over the 
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interval 1970-1971. This quiet time data set exhibits good uniform 
spatial distribution. The crustal field influnece is greatly attenuated 
at the spacecraft altitude of over 500 km and is assumed negligible for 
processing satellite data in the model. 

The observatory data was accomodated in the solution by solving for 
an independent constant anomaly bias magnetic field vector at the site of 
each observatory simultaneously with the spherical harmonic 
coefficients. A global set of 148 observatories was selected for the 
data set based on continuity of operation over the desired interval, data 
accuracy and spatial uniformity. The time series of annual means were 
examined for each observatory and jumps in the data caused by station 
relocation, equipment modification, etc. were identified and the 
intervals treated as independent stations with respect to the crustal 
biases. This resulted in a total set of 167 independent observatory 
vector biases from the 148 selected observatories. 

The repeat station data was selected based on data accuracy, spatial 
distribution and number of site occupations. At least three site 
occupations were required to qualify as acceptable. As the data set was 
generally sparse (and most often did not contain all three vector 
components), vector anomaly biases were not estimated for this data 
type. Instead, quadratic polynomial fits were made to the component time 
series and time derivative! (which are assumed to be independent of the 
crustal field) taken to be used as data for the main field model. 
Approximately 500 l3, (i, and t measurements from areas of Africa, South 
America and Australia were utilized. 

Marine survey data also contains contributions due to magnetic 
crustal anomalies. To obtain data for vast areas of the Atlantic, 

Pacific and Indian Oceans, forty-one long, straight tracks of surface 
marine scalar magnetic data (length greater than 1200 km) were selected 
and low pass filtered to remove wavelengths shorter than 500 km. This 
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filtered track of data was then sampled to pi’ovide marine data for he 
model data set. 

The weighting of the measurements in the least squares fit was based 
on the assumed accuracy of the data types. The satellite data was 
assigned a measurement noise sigma of 10 nT., while the sigma for 
individual magnetic observatories was determined from the scatter about a 
quadratic polynomial fit of the time series data. The invariant noise 
sigmas for the 167 observe cory vector components varied from 
approximately 5 to <10 nT. The sigma? used for the repeat station 
derivative data were also determined from polynomial fits, while tliat for 
the marine data (100 nT) was obtained by comparing filtered values of the 
long tracks at Intersecting pc' its. 

The PMAG (7/80) model has proven to be tlie best predictor of MAOSAT 
data among the available pre-MAGSAT models in the literature. The 
strength of the model is due to the extremely good POGO satellite 
coverage over the years 1965-1971 and the solution for observatory 
anomaly biases together with the inclusion of the third derivative 
spherical harmonic coefficients. It has been shown that the marine data 
and repeat station data has added no improvement to the predictive 
capability of the model. 

3.2 Mini -Batch Models 


The derivation of batch least squares main magnetic field models 
over approximatlely five year intervals from 1950-1976 was accomplished 
using the techniques described in Section 3.1 for PMAG (7/80) using the 
same POGO satellite data set and the same set of 167 magnetic 
observatories extending back to 1950. However, the filtered marine data 
and differentiated repeat station data from PMAG (7/80) was not used in 
the mini-batches. Instead, a selection of surface and airborne survey 
data was utilized to improve the data distribution within each mini-batch 
Interval. Mo attempt was made to model the crustal anomaly influence in 
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the survey data. The measurement noise sigma given this data in the 
estimation was set to reflect the crustal uncertainty. 

Because of differences in the distribution of the data and different 
data types, the processing of the five mini-batches was not identical. 
Table 3.1 shows the number of measurements for each data type, the 
weighting used for each data type and the maximum degree of the field 
coefficients and first derivatives for each of the mini-batches. The 
first satellite data was taken in 1965 (OGO-2) and thus the time 
intervals for the third and fourth mini-batches were adjusted so that 
this data was included in the third. Without the global satellite data, 
the distribution of the remaining data is sparse in some areas and thus 
the accurate recovery of all spherical harmonic coefficients and first 
derivatives is difficult. For this reason, the mini-batches from 1950 to 
1960 did not solve for all coefficients and derivatives. It is believed 
that the aliasing caused by truncation of the field model was small. 

The end time of the last batch was extended to include the latest 
observatory data available to us. The satellite data during this 
interval was taken in 1970-71 so that the solution is heavily weighted 
toward the beginning of the interval. 

The weighting sigmas used for the different data types were selected 
based rpon the fit of the data to the field model. Although this is not 
a rigorous procedure for determining the data accuracy, it is a 
reasonable approximation to it. The weighting sigmas for the observatory 
data were determined by the deviation of the data from a polynomial fit 
for each observatory. Thus the weighting varied from one observatory to 
the next. 

In each of these mini-batches, no corrections were made for the 
external field. The bias due to the external field should be accounted 
for as part of the main field while the high frequency variations of the 
external field will be treated as noise. 
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Table 3.1 Sunnary of Minf-Batch Field Models 


Time Span 

Number of 
observatories 

Number of Observations 

Maximum Degree and 
Order of Field Model 

Observatory 
< vector 
components) 

Survey/ 

airborne 

Satellite 

Total 

accepted 

Number 

edited 

constant 

linear 

1950-1955 

67 

849 

6776 

0 

7625 

153 

9 

6. 

1955-1960 

111 

1335 

12658 

0 

13993 

194 

10 • 

7 

1960-1966 

142 

2145 

58075 

3559 

63779 

468 

13 

13 

1966-1970 

148 

1578 

13545 

52016 

67139 

127 

13 

13 

1970-1976 

139 

2292 

9596 

15916 

27804 

189 

13 

13 


/ 


Data Type 

Data Weighting 

Data Editing Threshhold 

Observatory: X,Y,2 

5-45 nT* 

5000 nT 

Survey/Airborne: B, H, Z 

300 nT 

5000 nT 

0 

1.0“ 

5“ 

I 

0.5" 

5“ 

Satellite: B 

10 nT 

5000 nT 


*0 determined by polynomial fit to observatory data. 
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3.3 Computation of State Noise Covan ance and Mini-Batch Weighting 


Figures 3.1 are line printer plots of the field model coefficients 
for the five mini-batchs. There is one plot for all coefficients of the 
same degree. The code used to identify the coefficient order is: 


0 = g 
A = g 
B = h 
C = g 


1.0 

1.1 

1.1 

1.2 


Z 



Also shown on each plot is the range of the standard deviations for the 
coeficients (evaluated at the beginning of the interval) as computed from 
the least squares covariance matrix (inverse of information matrix). In 
general it may be stated that only the two or four coefficients of the 
highest order (for each degree) had large sigmas: most of the coeffi- 
cients had sigmas which were closer to the smaller of the two numbers. 

The number in parentheses is an "eyeball" average of the standard devia- 
tions where the sigmas of the high order coefficients were excluded. 

This is the number which was used to determine the scaling of the 
information matrix (to be explained shortly). 


These plots were made for two purposes: to determine whether or not 

the computed standard deviation was a realistic approximation of the 
coefficient accuracy and to aid in estimating the state noise variances. 
In examining these plots, we discovered that it was generally possible to 
fit a smooth curve through the five (or less) separate lines for each low 
degree coefficient. However, the deviation of coefficients from that 
curve was considerably larger than the computed uncertainty of the 
coefficient (i.e. computed sigma). This was disturbing because it 
suggests that either the sigmas which were used in weighting the data 
were wrong or that aliasing (because of the truncated field model or 





Figure 3.1 Plots of Mini-Batch Field Model Coefficients 

Degree 1 . 
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Figure 3.1 (continued), Denree 2 
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Figure 3.1 (continued). Degree 3 
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:l:1gure 3.1 (continued), Degree 4 
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Figure 3,1 (continued). Degree 7 
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other modeling errors) was significant. At this time, it is not known 
which explanation is more probable. The possibility of an error in the 
software has been investigated and it does not seem likely. 

In any case, it was obvious that the covariance matrix computed for 
each mini-batch had to be scaled to account for these modeling errors. 
Scaling the entire covariance matrix by a constant may not generate a 
matrix with the same structure as the true error covariance but it is 
certainly better than no scaling at all. The scaling factors for each 
mini -batch were determined by using the plots to estimate the true error 
sigma for all coefficients of the same degree. In other words, the 
deviation of the plots from a smooth curve (and from plots of other field 
models) was assumed to be a measure of the true accuracy. This number 
was divided by the average computed sigma (listed in parentheses) to 
compute a scale factor for each coefficient degree and then these scale 
factors for the different degrees were averaged to compute an overall 
scale factor for each mini-batch (Table 3.2). In general, the scale 
factors for coefficients of different degrees within the same mini-batch 
were quite consistent and thus there is some confidence in the method. 
However, this procedure for estimating the scale factors is subject to 
rather large uncertainty. 


Table 3.2 Estimated Scale Factors 



lini-Batci 


1950-55 

1955-60 

1960-66 

1966-70 

1970-76 

1.0 

4.0 

2.9 

2.0 

2.4 
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The second problem, determining the spectral density of the state 
noise, is somewhat more difficult to solve in a systematic manner. The 
deviation of the plots from a quadratic (or linear) curve should be a 
good indication of the state noise but since the uncertainty of the 
coefficients is so large, this deviation cannot be estimated accurately. 
Rather, a procedure which was based more upon the value of coefficients 
was used.- For the coefficients which were assumed to obey a quadratic 
law ( degrees 1 to 7), the standard deviation of the coefficient will grow 
as QgAt^/5 , Therefore, the expected deviation of the coefficient from 
the quadratic curve after a period of 20 years was chosen for coeffi- 
^cients of each degree and was computed from this. Likewise, the 

• O 

same procedure was used for coefficie nts which obey a linear trend where 
the coefficient uncertainty grows as-v/o^At ^ /3 . Table 3.3 lists the 

assumed coefficient uncertainty after 20 years for each degree and the 
resulting spectral density (Q^)* 
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Table 3.3 Estimated State Noise 


Coefficient 

Degree 

Estimated Coefficient 
Deviation cfter 20 years 

Time 

Model 

Spectral Density 
of White Noise Input 

1 

120nT 

quadratic 

.0224nT^/sec^ 

2 

80nT 

quadratic 

.OlOOnT^/sec® 

3 

60nT 

quadratic 

.0056nT^/sec^ 

4 

40nT 

quadratic 

.0025nT^/sec^ 

5 

30nT 

quadratic 

.0014nT^/sec^ 

6 

20nT 

quadratic 

.0006nT^/sec^ 

7 

14nT 

quadratic 

.0003nT^/sec^ 

8 

lOnT 

1 i near 

.0376nT^/sec^ 

9 

8nT 

linear 

.0240nT^/sec^ 

10 

6nT 

linear 

.0136nT^/sec^ 

11 

4nT 

linear 

.0060nT^/sec^ 

12 

2nT 

linear 

.0016nT^/sec^ 

13 

l,5nT 

linear 

.0008nT^/sec^ 
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3.4 Filtering Results 


Table 3.4 lists the statistics on the measurement residuals 
(difference between observed measurement and computed measurement) using 
MAGSAT data taken on March 15, 1980 and predictions from various field 
models. March 15 data was selected because this was a magnetically 
“quiet" day. The data was corrected for MAGSAT attitude errors and the 
effects of external magnetic fields and the vector data was limited to 
±45° in latitude. A brief description of the various field models is 
given below. 

AWC75 - Developed by Peddie and Fabiano using data from 1965 to 1975. 

Uses constant coefficients to degree 12 and linear time terms to 
degree 12. 

WMC80 - World Magnetic Chart 1980. Developed by Barker and Barracough 
using data from 1950 to 1977. Uses constant coefficients to 
degree 12 and linear time terms to degree 8. 

PMAG (7/80) - Developed by GSFC using data from 1960 to 1977. Uses cubic 
time coefficients to degree 4, quadratic time coefficients to 
degree 6, linear and constant terms to degree 13. 

Mini -Batch (1970-1977) - last mini -batch as described in section 3.1. 

GSFC (9/80) - Developed by GSFC combining MAGSAT Nov. 5 and 6 data with 
the PMAG (7/80) data set. Same form of time model as PMAG. 

GSFC (2/81) - Developed by GSFC using MAGSAT (Nov. 5 and 6, 1979 plus 
March 15, 1980), POGO (1965-1971) and observatory data over the 
period 1950 to 1977. Same form of time model as PMAG. 

Best filter model - field model from filter computer run which had 
highest likelihood function. 


45 






1 . 1 -f i .4 A' 







Table 3.4 Residual Statistics (nt) of March 15, 1980 MAGSAT Data for Various Field Models 

(corrected for attitude and external field) 


Geomagnetic Field Model 


Component 


AWC75 


B mean 67 
RMS 148 



PMAG(7/80) 

Mini-Batch 

1970-1976 

6SFC(9/80) 

GSFC(2/31) 

Best Filter 
Model 

(Pre-MAGSAT) 

-20 

5 

-3 

-1 

19 

92 

101 

11 

9 

78 

-9 

9 

-3 

-2 

11 

70 

no 

10 

8 

67 

2 

0 

1 

1 

2 

77 

119 

13 

12 

73 

9 

8 

-1 

-2 

-4 

109 

186 

13 

9 

108 

89 

127 

12 

9 

82 
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The mean and root-mean-square residual for the scalar and vector 
data are listed for each field model. The RMS value listed at the 
bottom is an RMS for all the measurements (approximately twice as many 
scalar measurements were used compared to vector measurements). 

Our filtered model is substantially better (8 to 39%) than any of 
the tested pre-MAGSAT models for prediction time intervals of 4 to 5 
years. (The date of the last data used in the WMC80 model is not 
precisely known). This demonstrates that the filtering method is clearly 
superior to the methods which have traditionally been used in the past to 
develop field models. 

The two GSFC models which used MAGSAT data were included in the 
table simply to show the limits of batch field modeling when little or no 
prediction is involved . (The GSFC (9/80) model was predicted 4 months). 
The final column lists the residual statistics for the filtered model 
which has been updated with MAGSAT data taken on November 5 and 6. 

Notice that this model is somewhat better than the GSFC (9/80) model and 
comparable to the GSFC (2/81) model (which used March 15 data). 

Table 3.5 lists the results of six filtering runs in which various 
inputs were perturbed to determine the optimum field model. Run #6 has 
the highest likelihood function and nearly the best fit to the data. The 
following list summarizes the results. 

1) The use of likelihood function in conjunction with the 
sum of weighted residuals was a reliable indicator of 
the field model which produced the best fit to the 
MAGSAT data. This is a significant result because it 
demonstrates that the filter model which is best 
supported by the data is also the best predictor of 
future data. 

(2) The best results were obtained when the nominal mini- 
batch scale factors were increased by 65%. See runs 
1, 2, 5 and 6. 
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Table 3.5 Summary of Field Models Produced by Kalman Filtering 


Residual Statistics (MA6SAT, March 15, 1980) 


Run #1 Scaling {State Noise 


Mean (nT) 


IxJPc^x 


i i jlikelihood 



note la 1 note 2a 


note la 


note 2b 


note la 


note 2c 


note la 


note 2d 


note lb I note 2a 


note lb 


note 2b 



(1) Scaling of mini-batches 

a) Nominal scaling (1.0, 4.0, 2.9, 2.0, 2.4) 

b) 1.65 X nominal scaling (1.7, 6.6, 4.7, 3.4, 4.0) 


(2) State Noise Spectral Density 
^S "" ^S,nom 

^ ^ ^S nom multiplied by 2) 

c) Qg = 16 (sigma multiplied by 4) 


d) = 64 Qg (sigma multiplied by 8) 


-1985 


lO 73 
C 
> 
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(3) The best results were obtained when the state noise 
spectral density listed in Table 3.2 was increased by 
a factor of 4 (the standard deviation multiplied by 
2). See runs 5 and 6. Attempts to vary the scaling 
of Qg between the linear and quadratic terms did 
not improve the results (early runs not listed). 

(4) The use of a mini-batch epoch time other than at the 
middle of the data span will degrade the results. 

This is supported by the results of runs not listed 
and runs made on simulated data. Presumably this is 
the result of aliasing errors (caused by the use of a 
linear time model in the mini-batch) which are a 
minimum at the center of the data span. 

In addition to the use of the log likelihood function to determine 
the optimum run, the sum of weighted residuals. 


can also be useful as an indicator of modeling errors. If all models are 
correct, this sum should be Chi-squared distributed with degrees of free- 
dom equal to the number of measureme.its, i.e. it should be approximately 
equal to the number of measurements. Since our measurements are actually 
the mini-batch field model coefficients, m should be equal to the sum of 
the number of coefficients for the five mini-batches, i.e. 


m = 147+183+390+390+390 = 1500. 
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In reality the sum of weighted residual should be slightly less than 1500 
because the a priori sigmas we used were considerably larger than the 
actual errors in the a priori estimate of the field model. It is diffi- 
cult to determine exactly what this sum should equal but a value of 1200 
seems reasonable. Notice that run #6 has the highest likelihood function 
and the sum of weighted residuals for this run is close to 1200. We 
found this sum to be quite useful in deciding how to scale the individual 
mini-batch covariance matrices and the state noise spectral density. 

The use of only the likelihood function to select the optimum model 
parameters can sometimes be misleading if the true measurement noise 
standard deviation is not known. This is demonstrated in Figure 3.2 
where the likelihood function and the RMS fit to the MAGSAT data is 
plotted for the six runs listed in Table 3.5. Using the nominal scaling 
of the mini-batches, the likelihood function is maximized when the state 
noise variances are multiplied by a number between 16 and 64. However, 
the sum of weighted residuals for these runs is much larger than 1200 and 
the prediction of the MAGSAT data is poor. When the scaling of the 
mini -batch variances is increased, a smaller scaling of the state noise 
variances will produce a higher value of the likelihood function (run 
6). Furthermore, the sum of residuals is closer to 1200 and the 
prediction of the MAGSAT data is improved. Because of the expense of 
making these computer runs, we did not attempt to find the optimum 
scaling parameters but run 6 is probably quite close to optimum. 

Table 3.6 lists the predicted set of field model coefficients at 
1980.0 for run 6. This model is called Geomagnetic Recursive Information 
Model 1977 (GRIM77). The degree and order listed in the table have been 
increased by 1 because Fortran will not allow zero subscripts. The field 
model listed in Table 3.7 (GRIM80) is similar but the results of run 6 
have been optimally combined with a field model based only upon November 
5 and 6, 1979 MAGSAT data. Thus Table 3.7 represents our best estimate 
of the geomagnetic field and its secular variation at epoch 1980.0. 
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J GRIM80 Field Model at Epoch 1980.0 (Includes MAGSAT Data) 
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4.0 DISCUSSION 


The results of the preliminary study show that recursive estimation 
has great potential for geomagnetic field modeling but more needs to be 
done to realize that potential. There are several problems in modeling 
that need to be addressed in order to achieve optimal estimation. 

First there is the question of the scaling of the covariance 
matrices produced in the mini-batches. When comparing the field 
estimates produced for different time spans, it is obvious that the 
actual errors -in the field coefficients are larger than the computed 
standard deviations. The cause of this modeling error is not known but 
it may be the result of incorrect weighting of the various measurement 
types. For example, it is possible that nonlinear weighting of the 
angular survey data as a function of the local field strength may be more 
appropriate. Other possiblities include aliasing caused by truncation of 
the spherical harmonic expansion {and secular variation terms) and 
program errors. These explanations and others must be examined until the 
cause of the discrepancy between the actual errors and the covariance 
matrix are resolved. It should not be necessary to scale the output of 
any mini-batches. This scaling is a crude attempt to account for 
modeling errors which should not be present. 

The determination of the optimum state noise variances is another 
problem which should be addressed in future work. Although we tried 
using different factors to scale the linear and quadratic state noise 
variances, the individual variances for each degree should be determined 
by a more systematic procedure than the eyeball estimate. In fact, it 
would be desireable to determine the variances for each degree and order 
separately, but this is probabaly an unrealistic goal (because of the 
computational burden). In our proposal, we suggested using maximum 
likelihood estimation (MLE) to determine the state noise variances. 
However, more study is required to find computationally efficient methods 
of performing MLE. 


PRECEDING PAGE BLANJT NOT FJLrviED 


BvsiyKss iM) TEcmoLoonuL Systems, Inc 


Other factors which should be examined in attempting to improve the 
field model include the use of marine and repeat survey data and the 
handling of observatory biases. It is questionable whether the use of 
marine and repeat survey data will improve the field prediction but it 
should be included in the mini-batches. A more significant factor is the 
handling of observatory biases. In the PMAG solution, one set of biases 
was estimated for the period 1960 to 1977 but a different set of biases 
was estimated for each mini-batch. Since the estimated biases were not 
entirely consistent from one mini-batch to the next (the estimated biases 
for the first two mini-batches were wildly inconsistent), the field 
estimates were likewise not consistent. An attempt was made to constrain 
the biases used in the mini-batches to the PMAG estimates but the 
estimate of the field model was further degraded and the covariance 
matrix was too small. Apparently there is some inconsistency between the 
PMAG bias estimates and the data in the individual mini-batches. It 
would be preferable to combine the mini-batch bias estimates using the 
filter to produce an estimate of the biases for the entire period. 
Unfortunately this would greatly increase the computational burden, both 
in the mini-batches and in the filtering. 

Finally, the lack of consistency in the field coefficient estimates 
for the 1955 to 1960 mini-batch compared to the other mini-batches may 
have degraded the final filtered estimate. The inconsistency of this 
mini-batch may have resulted from over-parameterization of the field 
model (use of a maximum coefficient degree which is too high) for the 
available data. Thus minor modeling errors such as incorrect weighting 
of the measurements could have caused highly erroneous bias and field 
model estimates. This explanation can be tested easily and should be 
investigated in future studies. 
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5.0 CONCLUSIONS AND RECOMMENDATIONS 


5,1 Conclusions 


(1) The feasibility of using Kalman filtering for geomag- 
netic field modeling has been clearly demonstrated. 

The prediction of the MAGSAT data using the filtered 
field model was significantly better than the predic- 
tion of other tested pre-MAGSAT field models. This 
was accomplished even though some modeling errors 
existed in the mini -batch field models. 

(2) The use of the likelihood function in conjunction with 
the sum of weighted residuals to determine the optimum 
models has been validated. There was a strong 
positive correlation between the value of the 
likelihood function and the ability of the field model 
to predict the MAGSAT data. 

(3) It is believed that the dominant error source in 
performing the mini-batch data reduction was incorrect 
weighting of the survey measurements. 


5.2 Recommendations 


(1) The source of the modeling errors in the mini-batch 
data reduction should be determined and eliminated. 
This is a necessary step for further filtering work 
and for any type of geomagnetic field modeling. 

(2) Marine and repeat survey data should be included in 
the mini -batch data reduction. 
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(3) It would be desireable to solve for a common set of 
observatory biases for all the mini-batches using the 
Kalman filter. This should eliminate some of the 
inconsistency between field estimates for different 
mini-batches. 

(4) After the modeling errors have been eliminated, 
optimum values of the state noise variances should be 
determined. Ideally this should be performed using 
maximum likelihood estimation and may require the use 
of data back to 1900. 

(5) Field models for times in the past should be obtained 
using an optimal smoother based upon the results of 
the filter. 

(6) Contour maps showing the accuracy of the field model 
can be produced using the error covariance matrix 
computed by the filter. 
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Appendix A.l 

COMPARISON OF KALMAN AND 
WIENER FILTERING 


The review of our proposal for this study requested that the Kalman 
filtering approach be compared to Wiener filtering. In general, it may 
be stated that Wiener filtering is a subset of Kalman filtering: Wiener 

filtering requires certain assumptions that are not required in Kalman 
filtering. It may appear that Kalman filtering requires knowledge of the 
system dynamics but, in fact, no more knowledge is required than for Wiener 
filtering [1]. 

In this section, an explanation of Wiener filtering, a discussion 
of how the geomagnetic field problem fits into the framework of Wiener 
filtering, and a comparison between Wiener and Kalman filtering for this 
specific problem will be given. Other references which discuss, more generally, 
the comparisons between the two approaches are [l ,4 ]• 

The basic problem addressed in Kalman and Wiener filtering is one 
where we have observations z(t) of a signal process y(t) corrupted 
by additive white noise v(t) 

z(t) '■ y(t) + v(t) tQ < t < t^ (!) 


where 


E(v(t)v^(T)) = V 5(t-x) . 

Although Wiener and Kalman filtering can handle the case where y(0 is 
correlated with v(*) , (e.g. this occurs when feedback is used) we shall 
only consider the uncorrelated case 

E(y(t)v'^(T)) 5 0 . 
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It should be noted that Wiener's original work was limited to the case 
where z , y and v were scalars, y(t) and v(t) were statistically 
stationary and the observation interval was infinite, (t^ = -»). Although 
his work has been extended to avoid these restrictions, the following descrip 
tion of Wiener filtering is limited to scalar observations over an infinite 
interval . 



y{-r)=y(t) 


Figure 1 General Filtering/Prediction Problem 

Figure 1 shows the general filtering/ smoothing/prediction problem. 
We wish to find a function h(t,t} in 


y(t|t^) = h{t,T)z{T)dT 


such that tr E[y(t)-y(t)][y(t)-y(t)] is minimized. The three separate 
problems are defined as; 


filtering : 
smoothing : 
prediction: 


t = t^ ■ 
to < t < t^ 
t > tj- . 


Although we are primarily interested in prediction of the geomagnetic field, 
we shall confine the present discussion to filtering (t=t^) since the 
derivation for prediction is similar. 


It is well known [ 3] that the least squares solution to (2) is obtained 
by the orthogonality property 
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E[{y(t)-y{t))z”'’(T)] = 0 tg < T < t . (3) 

Using (2) and (3) we find that 
t 

h(t,T)V + I h(t,x)E[y(^)y"^(T)]dx = E[y(t)y^(i:)] Cq - ^ 

We now apply Wiener's assumptions that tg = -» and E[y(t)y^(T}] is statis- 
tically stationary, i'.e. depends only on jt-T| . Then by making a change 
of variables, equation 4 becomes 

h(t)V + h{x) ^y(t-X)dX = <i.^(t) 0 < t < » ( 5 ) 

where = E[y(t)y^(0)] = E[y(-t)y^(0)] . Equation (5) is the famous 

Wiener-Hopf equation. Wiener solved this using spectral factorization 
(for the scalar case). For the processes with rational spectral densities, 
^y(t) has the form 

♦y(t) = fy(ltl) = a, (6) 

where a^. and are constants (possibly complex). The bilateral Laplace 
transform of (6) 

<^y(s) = <}>y(t)e"^^dt (7) 

2 

is a ratio of polynomials in s whose poles always occur in complex 
conjugate pairs: every root in the right-half s-plane is accompanied 

by a root in the left-half plane. Thus, 4)(s) can be written as 

%(s) = «y(s) ‘&y(-s) (8) 

where 4'y is a unique factor of $(s) . Using this factorization, Wiener 
showed that the convolution in equation (5) could be written in terms of 
these spectral factors and the optimum H(s) is found to be 
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H(s) 


V+$y(s) 


( 9 ) 


A brief description of the Kalman filter was given in preceding sections 
of this report. We shall not elaborate further on the details of Kalman 
filtering since they are generally well known and numerous fine references 
exist on the subject [for example, see 2 or 3 ]. 

In comparing the development of the Wiener filter with the Kalman 
filter, several differences in assumptions are apparent. 

(1) The (original) Wiener filter assumes stationarity of 
y{t) and v(t) while the Kalman filter makes no such 
assumption. As we have noted, Wiener filtering has been 
extended to handle non-stationarity but the procedures 
tend to be complicated. 

(2) The (original) Wiener filter assumed that the signal 
z(t) was available from time = while the Kalman 
filter operates on finite data spans. The Kalman filter 
can be made equivalent to a Wiener filter by using the 
steady state Kalman gains but the Wiener filter is not 
so easily modified to handle finite data spans. This 
limitation of the Wiener filter is a problem for geo- 
magnetic modeling because the data span is relatively 
short compared to time constants of the field dynamics. 

(3) It would appear that the Kalman filter may only be used 
when the signal process is the output of a known finite- 
dimensional system driven by process noise. In contrast, 
the Wiener filter makes no such assumption. Instead, it 
requires knowledge of the covariance functions E(y(t)y (t)] 

T 

and E(z(t)z(x)] . However, [1 ] shows that a Kalman 
filter may be constructed using only the same covariance 
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functions used in the Weiner filter (although the 
state error covariances computed by the Kalman filter 
may not have any physical meaning). Furthermore, the 
solution of the Wiener-Hopf equation requires that the 
system be finite-dimensional. Thus, the distinction 
between the two approaches is not as significant as 
might be presumed. 

The implications of these different assumptions for geomagnetic field 
modeling are quite significant. Particularly for this problem, the Kalman 
filter is to be preferred to the Wiener filter. The stationarity assump- 
tion used in the Wiener filter is certainly not valid for geomagnetic 
modeling. The "measurements" input to the filter must be the mini -batch 
estimates of the. geomagnetic field rather than the raw magnetic field 
measurements obtained at specific points on the earth. (If the thousands 
of raw measurements were used directly in the filter, the computational 
burden would be enormous). Thus the covariance of v(t) will vary from 
one mini -batch to the next. Furthermore, the signal process is not stationary 
for the limited time period of our data. Since the data span is short 
compared to the time constants of the geomagnetic field dynamics, we have 
modeled (in the Kalman filter) the dynamics as free (no feedback) integrators 
fed by white noise. Obviously the output of these integrators will not 
be bounded or stationary and thus the model is only approximately valid 
for relatively short periods of time. This is a reasonable approximation 
for a Kalman filter, but the lack of stationarity is a serious problem 
for the Wiener filter. 

The extremely short data span (5 time points) is also a serious draw- 
back to the use of the Wiener filter. This has two bad effects; the 
sample covariance function E[z(t)z^(x)] will be grossly in error (and 
truncated) and the optimal weighting function h(t) will not have a sufficient 
amount of data to operate on. The only assumption used in the Kalman 
filter which may cause problems is the finite dimensional dynamic model. 

If this model vjere grossly in error, then the estimates produced by the 
filter would be erroneous. However, this is not a serious problem for 
geomagnetic modeling since we are sampling a short time span of a system 
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with long time constants. As we have already pointed out, the alternative 
(use of a sample covariance function) is much lass’ attractive. 

Although it is possible to eliminate some of the restrictive assump- 
tions imposed upon Wiener filtering, the methods used to overcome the 
restrictions are not simple (particularly for the non-stationary, multi- 
dimensional, discrete time, limited data span case). In contrast, Kalman 
filtering seems ideally suited to the task. We are aware of no advantages 
(either analytically or computationally) that Wiener filtering enjoys 
relative to Kalman filtering for this application. 
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Appendix A. 2 

Discussion of the Predictive Properties of 
Kalman Filtering Versus Batch Least Squares 

In a review of a proposal for geomagnetic modeling by Kalman filtering, 
one reviewer made the statement that "Kalman filtering is known to be not 
so good for prediction." Although there are some problems for which the 
long term predictions of Kalman filters are inferior to those of batch 
processors, the preceeding statment is generally not true . As evidence 
of this, numerous applications papers exist demonstrating (successfully) 
the predictive ability of Kalman filters. It is the purpose of this appendix 
to discuss the conditions under which the Kalman predictor is inferior 
or superior. 

Consider the case where measurements of a dynamic process are to be 
input to a Kalman filter and a batch least squares estimator where the 
same measurement and dynamic models are used in both estimators. The only 
difference between the two estimators is the inclusion of process (state) 
noise statistics in the Kalman filter. We also assume that the dynamic 
model used in the estimators is not a perfect match to the real world model. 
Thus both estimators will produce state estimates which contain errors 
due to the modeling errors. We restrict our attention to the state estimate 
produced at the end of the data span so that both estimators have processed 
an equal amount of data. Under the assumption that the state noise covariance 
used in the Kalman filter is a reasonable approximation to the actual modeling 
errors in the dynamics, the Kalman filter will almost always produce a 
better estimate of the final state (at the end of the data span) than the 
batch estimator. In fact, the estimation errors of the batch estimator 
tend to be largest at the ends of the interval . Thus, the Kalman predictions 
based upon the final state estimate will usually be better than those of 
the batch estimator for "short term" predictions. However, under some 
circumstances, the long term predictions of the batch estimator will be 
better than those of the Kalman filter. This is particularly true when 
the modeling errors tend to be cyclic. 
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As an example of this, consider satellite orbit determination. As 
with most real world problems, it is almost impossible to model the system 
perfectly. Typical orbital modeling errors include solar radiation, atmos- 
pheric drag, gravity field errors, tracking system biases and refraction 
effects. When orbit predictions from Kalman filters are compared with 
those of batch estimators, it is usually found that the error in the Kalman 
estimate is greater than that of the batch estimator for prediction intervals 
greater than one quarter orbit. For prediction less than one-quarter orbit, 
the Kalman filter is usually superior. This can be explained on the basis 
of the batch processor's global fit to the data: it is attempting to average 

the modeling errors over the entire data span, whereas, the Kalman filter 
applies the most weight to data at the end of the interval. In other words, 
the batch processor behaves somewhat like a low pass filter and treats the 
dynamic errors as high frequency noise with zero mean . Since the propogation 
of orbital errors is cyclic (because of orbital dynamics), the prediction of 
the batch processor usually tends to match the long term behavior of the orbit. 

In order to better understand this phenomenon, consider a trivial 
example which contains large modeling errors. Assume that we have measure- 
ments of z(t) = oostii^t (0 < t < — ) and that we model the process as a 

a - - (Uq 

constant plus linear term: 

z = X-| -f V 



X2 = q 


where v and q are white noise processes. We use this same model in 
both the Kalman filter and batch least squares estimator (although q is 
obviously assumed to be zero in the batch estimator). We also assume that 
the contribution of a priori information is negligible so that the state 
estimates are dominated by the measurement data. It is aasy to show that 
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the estimates of x-j and X 2 from the batch estimator will be identically 
equal to zero. However, the estimates from the Kalman filter depend upon 
the value of Q^5(t-T) = E(q(t)q (t)) used. Figure 1 shows typical behavior 
of the Kalman filter. 



Figure 1 Input/Output of Kalman Filter and Batch 
Least Squares Estimator for Example 1 


Notice that for medium values of , the Kalman filter estimate 


of x-j at 


t = 

“0 


is fairly accurate but the prediction rapidly diverges 


from the true z(t) . For higher values of , both the estimate and 
the derivative are more accurate. For low values of , the estimate 
of x-| and X 2 are closer to those of the batch estimate (i.e. zero). 
Notice the batch least squares estimator produces better long tern predic 
tions than the Kalman filter. 


Now consider a minor modification of this example which is closer 
to the problem of geomagnetic modeling. Instead of starting at time 
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t = 0 , start at t = tt/uiq . Then we get the results shown in Figure 2 



Figure 2 Input/Output of Kalman Filter and Batch 
Least Squares Estimator for Example 2 

Here the prediction of the batch estimator is completely erroneous 
while the Kalman prediction is much closer to the actual z(t) . However, 
even the Kalman filter cannot accurately predict the signal because of 
the large modeling errors. 

These examples were intentionally chosen to exaggerate the best and worst 
properties of the Kalman and least squares estimators but, in many respects, 
these examples are typical of real world problems. As we have noted, the first 
example is similar to the orbit determination problem (usually one revolution 
or more of data will be processed) while the second example is somewhat 
similar to geomagnetic field modeling. The second example is also more 
typical of problems where the underlying process really is stochastic (or 
appears to be stochastic) rather than deterministic. When attempts are made to 
predict the behavior of stochastic signals, then it is important for the 
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estimator to have the best possible state estimate at the ^ of the data 
span. That is exactly what the Kalman filter is intended to do. 
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ADoendix A-3 Plots of WC 80 Field Model Coefficients 


The following line printer plots of the WC80 field model were 
produced using almost the same scale as the plots of the mini-batch 
coefficients (Figure 3.1). Directly above each plot are four numbers. 

The first number is the vcoefficient degree. The second and third numbers 
are the minimum and maximum values of the plot in nanotesla. The last 
number is the scale factor in nT per small division (one print position). 

The coding of the coefficient order is the same as in Figure 3.1, 


The numbers at the top of each plot are the print positions for each 
letter (order) for the ten points plotted. This is useful in resolving 
plotting ambiguities. 
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The following line printer plots of the PMAG (7/80) field model were 
produced using almost the same scale as the plots of the mini-batch 
coefficients (Figure 3.1). Directly above each plot are four numbers. 

The first number is the coefficient degree. The second and third numbers 
are the minimum and maximum values of the plot in nanotesla. The last 
number is the scale factor in nT per small division (one print position). 

The coding of the coefficient order is the same as in Figure 3.1, 


The numbers at the top of each plot are the print positions for each 
letter (order) for the ten points plotted. This is useful in resolving 
plotting ambiguities. 
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Appendix A-5 Plots of 6SFC (2/81) Field Model Coefficients 

The following line printer plots of the GSFC (2/81) field model were 
produced using almost the same scale as the plots of the mini-batch 
coefficients (Figure 3.1). Directly above each plot are four numbers. 

The first number is the coefficient degree. The second and third numbers 
are the minimum and maximum values of the plot in nanotesla. The last 
number is the scale factor in nT per small division (one print position). 

The coding of the coefficient order is the same as in Figure 3.1, 

i .e. 


0 

A 

B 

C 


= h'-i 


The numbers at the top of each plot are the print positions for each 
letter (order) for the ten points plotted. This is useful in resolving 
plotting ambiguities. 




ASX 


43 CS. 


iS)C3. 


0 

1 



u_ 

>0 

Vi 

■> 

o 


o 

o 


ORieiNAU PAGE B 
OF POOR QUALITY 


rvi 


r- 

r* 

tn 


in 

o 


ui 

^iDknaiinuiininaiMiin o 




















taasasBf 


Q 

A 

t 

C 

0 

E 

F 

G 

II 

1 

J 

K 

L 

8j 

79 

6E 



19 

72 

69 

66 

69 

67 

52 

65 

D1 

76 

67 

69 

aa 

20 

76 

69 

64 

68 

67 

51 

65 

at 

71 

67 

69 

BB 

20 

76 

69 

64 

66 

67 

51 

65 

79 

70 

67 

69 

39 

22 

79 

69 

62 

68 

66 

49 

66 

79 

7E 

67 

69 

B9 

22 

79 

69 

62 

68 

66 

49 

66 

78 

£G 

67 

70 

B9 

23 

82 

69 

61 

68 

66 

49 

67 

70 

80 

67 

70 

U9 

23 

B2 

69 

61 

68 

66 

49 

67 

78 

at 

66 

72 

89 

26 

83 

69 

00 

69 

66 

48 

68 

78 

B1 

66 

72 

89 

26 

83 

69 

60 

69 

66 

48 

68 

78 

81 

66 

79 

8B 

28 

83 

69 

60 

70 

67 

47 

70 

t 

- 

339.9 


263.4 


5 

.0 27 



















0 

93 

A 

56 

C 

10 

6? 

l2 

sl 

G 

56 

H 

69 

I 

67 

J 

80 

K 

76 

L 

60 

n 

72 

N 

61 

S3 

54 

64 

10 

6} 

73 

6B 

58 

70 

60 

19 

75 

60 

72 

62 

S3 

54 

54 

10 

61 

73 

60 

58 

70 

68 

79 

75 

60 

72 

62 

93 

53 

51 

10 

60 

73 

68 

59 

71 

611 

79 

75 

60 

71 

63 

91 

53 

51 

10 

60 

73 

6B 

59 

71 

63 

79 

75 

60 

71 

63 

S3 

52 

4S 

10 

60 

74 

68 

61 

72 

69 

78 

75 

61 

71 

64 

9.3 

52 

49 

10 

60 

74 

68 

61 

72 

69 

78 

75 

61 

71 

64 

94 

52 

46 

10 

60 

75 

68 

62 

73 

69 

77 

74 

61 

70 

64 

94 

52 

46 

10 

60 

75 

68 

62 

73 

69 

77 

74 

61 

70 

64 

94 

50 

43 

10 

60 

76 

68 

64 

74 

70 

76 

74 

62 

69 

66 


- 202.2 


1b1.lt 


2 . * 11 )? 





0 

A 

6 

C 

D 

t 

J 

6 

K 

L 

J 

K 

L 

H 

H 

n 

l> 

(i2 

64 

65 

60 

59 

|9 

64 

63 

59 

66 

64 

58 

70 

69 

65 

61 

57 

63 

64 

65 

61 

59 

5S 

64 

62 

55 

66 

64 

58 

70 

68 

64 

£1 

57 

63 

64 

65 

61 

59 

59 

64 

62 

59 

66 

64 

58 

70 

68 

64 

61 

57 

6U 

64 

6 5 

61 

56 

59 

64 

62 

5fl 

65 

64 

59 

69 

67 

63 

61 

57 

64 

64 

65 

61 

58 

55 

64 

62 

58 

65 

64 

59 

69 

67 

63 

61 

57 

(5 

64 

64 

61 

sa 

59 

64 

61 

57 

65 

64 

60 

68 

66 

62 

£1 

58 

65 

64 

64 

61 

sa 

55 

64 

61 

57 

65 

64 

60 

68 

66 

62 

61 

58 

f>6 

64 

64 

62 

58 

59 

63 

61 

57 

64 

64 

61 

68 

65 

61 

62 

58 

66 

64 

64 

62 

50 

59 

63 

61 

57 

E4 

64 

61 

68 

65 

61 

62 

58 

67 

64 

64 

62 

57 

59 

63 

60 

56 

63 

64 

62 

67 

65 

59 

£2 

58 


8 -2U.9 


209.7 


3,555 










;; 

*2 

A 
t 1 

0 

48 

C 

£9 

6? 

54 

P 

60 

R 

63 

II 

59 

I 

61 

J 

59 

K 

sa 

ol 

B 

58 

II 

64 

U 

59 

P 

6 1 

bS 

n 

56 

s 

64 

62 

48 

£9 

65 

54 

60 

63 

59 

60 

58 

58 

63 

59 

64 

59 

61 

62 

57 

? 

e>i 

02 

4 6 

£9 

OS 

54 

60 

63 

59 

60 

58 

58 

63 

59 

64 

£9 

6i 

62 

57 

§ 

£4 

62 

49 

£9 

65 

54 

61 

63 

58 

00 

58 

58 

61 

59 

64 

59 

60 

61 

57 


£4 

02 

49 

£9 

65 

54 

61 

63 

58 

60 

58 

58 

63 

59 

64 

59 

60 

61 

57 

'■j 

63 

62 

49 

£9 

6 5 

54 

61 

63 

5E 

£9 

57 

58 

63 

60 

64 

59 

59 

6 0 

58 


f1 

62 

49 

£9 

65 

54 

61 

63 

5£ 

59 

57 

58 

63 

60 

64 

£9 

59 

60 

58 

t 

bj 

03 

45 

£9 

65 

54 

62 

63 

56 

ca 

£7 

58 

63 

60 

64 

59 

50 

59 

59 

1 

63 

63 

49 

£9 

65 


62 

61 

58 

59 

57 

58 

63 

00 

64 

59 

58 

59 

59 

< 

62 

6 3 

50 

£9 

66 

£4 

62 

63 

57 

58 

56 

58 

63 

61 

fi 

60 

57 

58 

59 


S -13T.1 146„7 2.165 





















V7'.5JE»5>'’ 

u'C>4'>;tv ■ 





n 

K 

G 

r 

I) 

B 

r 

8 

!l 

T 

.1 

K 

L 

B 

H 

0 

P 

Q 

R 

s 

T 

0 

V 

7fl 

6« 

bll 

61 

RO 

91 

68 

‘F 

59 

63 

61 

65 

61 

70 

78 

75 

62 

59 

4 7 

90 

72 

58 

111 

71 

b4 

69 

61 

79 

80 

67 

60 

59 

64 

69 

65 

6? 

71 

74 

7 b 

62 

59 

50 

HR 

70 

f.1 

104 

71 

fi*) 

69 

61 

79 

88 

67 

60 

59 

64 

69 

65 

62 

71 

74 

75 

62 

59 

50 

RR 

■fO 

61 

104 

77 

(>4 

69 

61 

77 

85 

6S 

6? 

58 

64 

69 

65 

61 

71 

70 

75 

63 

60 

53 

R5 

fill 

65 

97 

77 

64 

69 

61 

77 

BS 

65 

62 

58 

64 

69 

65 

63 

71 

70 

75 

63 

60 

53 

H5 

eii 

65 

97 

7U 

64 

70 

61 

76 

82 

64 

64 

58 

65 

69 

66 

64 

72 

66 

15 

63 

61 

56 

R2 

6/ 

61 

“0 

7« 

64 

70 

61 

76 

82 

64 

64 

58 

65 

69 

66 

64 

72 

66 

75 

63 

61 

56 

H2 

67 

69 

90 

7S 

6 4 

71 

60 

7S 

79 

6? 

65 

58 

F5 

61 

66 

65 

72 

62 

75 

64 

62 

59 

79 

65 

72 

R3 

7S 

64 

71 

60 

75 

79 

62 

65 

58 

65 

69 

66 

65 

72 

62 

75 

64 

62 

59 

79 

65 

72 

R3 

77 

6J 

71 

60 

74 

72 

60 

68 

57 . 

66 

69 

66 

66 

7J 

58 

75 

65 

63 

63 

76 

6 3 

77 

74 

1 1 

- 

19.55 



16.10 


0. 

2171 




















c 

i|7 

it 

50 

.1 

,C 

<9 

s5 

5^ 

6l 

9§ 

a!l 

s5 

.1 

98 

58 

5^f 

9§ 

9§ 

0 

33 

5/ 

6? 

R 

40 

S 

62 

T 

59 

It 

59 

V 

55 

» 

16 

9^ 


9^ 

*19 

50 

97 

10 

58 

53 

67 

97 

90 

59 

9fi 

57 

56 

SO 

98 

37 

56 

69 

94 

59 

57 

55 

54 

23 

94 

22 

91 

It'} 

50 

97 

70 

SB 

S3 

67 

97 

90 

59 

98 

57 

56 

50 

90 

37 

56 

69 

99 

59 

57 

55 

59 

23 

99 

22 

91 

50 

50 

«e 

tb 

57 

52 

65 

99 

9'l 

59 

99 

55 

55 

52 

51 

40 

55 

61 

98 

57 

55 

55 

54 

29 

97 

29 

93 

50 

50 

98 

Eb 

57 


65 

99 

99 

59 

99 

55 

55 

52 

51 

90 

55 

61 

98 

57 

55 

55 


29 

97 

29 

93 

52 

5C 

5C 

62 

57 

51 

69 

51 

97 

59 

50 

54 

59 

53 

53 

99 

55 

59 

52 

55 

59 

55 

53 

35 

99 

36 

45 

52 

50 

50 

62 

5? 

£1 

69 

51 

47 

59 

50 

59 

59 

S3 

53 

99 

55 

59 

52 

55 

59 

ei; 

.53 

35 

49 

36 

95 

53 

50 

52 

58 

56 

50 

62 

53 

50 

50 

51 

53 

52 

59 

56 

98 

59 

56 

55 

52 

51 

56 

53 

41 

52 

93 

96 

53 

SO 

52 

58 

56 

50 

62 

53 

50 

50 

51 

53 

52 

59 

56 

90 

59 

56 

55 

52 

53 

56 

51 

41 

52 

43 

96 

55 

5C 

59 

59 

55 

99 

59 

56 

59 

. 50 

52 

51 

51 

55 

59 

52 

53 

53 

60 

99 

52 

56 

52 

48 

55 

51 

98 


13 -9.015 11.71 0.1727 



ORIGINAL PAGE !S 




